home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / demo2.i < prev    next >
Text File  |  1995-09-26  |  7KB  |  255 lines
  1. /*
  2.    DEMO2.I
  3.    Mesh plotting demo
  4.    To run, start yorick, then type the following two lines:
  5.  
  6.      #include "demo2.i"
  7.      demo2
  8.  
  9.    You can rerun the movies as many times as you want by typing "demo2".
  10.    To run only the first, second, or third movie, type "demo2, n"; for
  11.    example to run only the third movie, type
  12.  
  13.      demo2, 3
  14.  
  15.  
  16.    $Id$
  17.  */
  18. /*    Copyright (c) 1995.  The Regents of the University of California.
  19.                     All rights reserved.  */
  20.  
  21. func demo2(which, time_limit)
  22. /* DOCUMENT demo2
  23.      Exhibit quadrilateral mesh plots in 3 movies of a drumhead.
  24.      The drumhead is initially stationary, but has a bump near one
  25.      edge.  Yorick is solving a 2D wave equation to compute the
  26.      evolution of this bump.
  27.  
  28.      The first movie is a filled mesh plot with color "proportional"
  29.      to the height of the surface of the drum.  A few well chosen
  30.      contour levels (here 3) add a lot to a filled mesh plot.
  31.  
  32.      The second movie is a "3D" perspective plot of the height of the
  33.      drumhead.  In this movie, the mesh lines are drawn, which is
  34.      slightly confusing since the cells are not all the same shape.
  35.  
  36.      The second movie is a "3D" shaded plot of the height of the
  37.      drumhead.  Yorick computes surface shading based on the angle
  38.      of each cell from a light source.
  39.  
  40.      As you watch this, you might reflect on the two dimensionality
  41.      of your retina.  What Yorick lacks by way of 3D graphics is
  42.      really just fancy hidden surface algorithms; the simple
  43.      painter's algorithm used here and in plwf.i is easy to
  44.      implement.
  45.  
  46.      There are two optional arguments to demo2: the first is the
  47.      number fo the movie (1, 2, or 3) you want to watch; the second
  48.      is a time limit on the duration of each movie in seconds (default
  49.      is 60 seconds each).
  50.  */
  51. {
  52.   require, "movie.i";
  53.  
  54.   /* generate a 30-by-30 cell mesh on the [-1,1] square */
  55.   x= span(-1, 1, 31)(,-:1:31);
  56.   y= transpose(x);
  57.   /* map the square mesh into a mesh on the unit circle
  58.      this mesh has more nearly equal area cells than a polar
  59.      coordinate circle */
  60.   scale= max(abs(y),abs(x))/(abs(y,x)+1.e-30);
  61.   /* note that abs(y,x)=sqrt(x^2+y^2) */
  62.   x*= scale;
  63.   y*= scale;
  64.  
  65.   f= f0= exp(-8.*abs(y+.67,x+.25)^2)*(1.-abs(y,x)^2);
  66.   fdot= 0.0*f(2:-1,2:-1);
  67.  
  68.   af= abs(f(2:-1,2:-1));
  69.   lf= laplacian(f, y,x);
  70.   dt= 0.08*sqrt(2.*max(af)/max(abs(lf)));
  71.  
  72.   window, 0, wait=1, style="nobox.gs";
  73.   palette, "heat.gp";
  74.   limits, -1, 1, -1, 1;
  75.  
  76.   fma;
  77.   plt, "Yorick         (c) 1994 Regents of the U.C.\n"+
  78.     "X Division, Lawrence Livermore Lab\n"+
  79.       "\nInterpreter for scientific computing\n"+
  80.     "UNIX and MacIntosh versions\n"+
  81.       "anonymous FTP from:\n"+
  82.         "ftp-icf.llnl.gov:/pub/Yorick\n"+
  83.           "\nThis demo shows three movies\n"+
  84.         "of a drumhead oscillating",
  85.         0.40, 0.64, justify="CH", font="helvetica", height=24;
  86.   redraw;
  87.   pause, 8000;
  88.   fma;
  89.  
  90.   /* roll the filled mesh movie */
  91.   if (which && which!=1) goto persp;
  92.   fc= f(zcen,zcen);
  93.   cmin= cmax= max(abs(fc));
  94.   cmin= -cmin;
  95.   level= cmax/4.;
  96.   movie, display_plf, time_limit;
  97.   pause, 1000;
  98.   animate, 0;
  99.   write,format="%f frames of filled mesh drumhead completed in %f sec\n", 
  100.         movie_timing(4), movie_timing(3);
  101.   write,format="Rate for filled mesh is %f frames/(CPU sec), %f frames(wall sec)\n",
  102.     movie_timing(4)/(movie_timing(1)-movie_timing(5)+1.0e-4), 
  103.     movie_timing(4)/(movie_timing(3)-movie_timing(5)+1.0e-4);
  104.  
  105.   /* roll the perspective movie */
  106.   persp: if (which && which!=2) goto shade;
  107.   f= f0;
  108.   limits, -1, 1, -1, 1;
  109.   movie, display_plm, time_limit;
  110.   write,format="%f frames of wireframe surface drumhead completed in %f sec\n", 
  111.         movie_timing(4), movie_timing(3);
  112.   write,format="Rate for filled mesh is %f frames/(CPU sec), %f frames(wall sec)\n",
  113.     movie_timing(4)/(movie_timing(1)-movie_timing(5)+1.0e-4), 
  114.     movie_timing(4)/(movie_timing(3)-movie_timing(5)+1.0e-4);
  115.  
  116.   /* roll the shaded movie */
  117.   shade: if (which && which!=3) return;
  118.   f= f0;
  119.   limits, -1, 1, -1, 1;
  120.   movie, display_pl3, time_limit;
  121.   write,format="%f frames of filled surface drumhead completed in %f sec\n", 
  122.         movie_timing(4), movie_timing(3);
  123.   write,format="Rate for filled mesh is %f frames/(CPU sec), %f frames(wall sec)\n",
  124.     movie_timing(4)/(movie_timing(1)-movie_timing(5)+1.0e-4), 
  125.     movie_timing(4)/(movie_timing(3)-movie_timing(5)+1.0e-4);
  126.  
  127.   fma;
  128.   limits;
  129. }
  130.  
  131. func display_plf(i)
  132. {
  133.   /* display first */
  134.   fc= f(zcen,zcen);
  135.   cmin= cmax= max(abs(fc));
  136.   cmin= -cmin;
  137.   plf, fc, -y, -x, cmin=cmin, cmax=cmax;
  138.   /* the 0 contour level is too noisy without some smoothing... */
  139.   fs= f(zcen,zcen)(pcen,pcen);
  140.   plc, fs, -y, -x, levs=0., marks=0, color="green", type="solid";
  141.   plc, f, -y, -x, levs=level, marks=0, color="black", type="dash";
  142.   plc, f, -y, -x, levs=-level, marks=0, color="green", type="dash";
  143.  
  144.   /* then take a step forward in time */
  145.   lf= laplacian(f, y,x);
  146.   af= abs(f(2:-1,2:-1));
  147.   fdot+= lf*dt;
  148.   f(2:-1,2:-1)+= fdot*dt;
  149.  
  150.   return i<200;
  151. }
  152.  
  153. func display_plm(i)
  154. {
  155.   /* display first */
  156.   pl3d,0, f, y, x;
  157.  
  158.   /* then take a step forward in time */
  159.   lf= laplacian(f, y,x);
  160.   af= abs(f(2:-1,2:-1));
  161.   fdot+= lf*dt;
  162.   f(2:-1,2:-1)+= fdot*dt;
  163.  
  164.   return i<200;
  165. }
  166.  
  167. func display_pl3(i)
  168. {
  169.   /* display first */
  170.   pl3d,1, f, y, x;
  171.  
  172.   /* then take a step forward in time */
  173.   lf= laplacian(f, y,x);
  174.   af= abs(f(2:-1,2:-1));
  175.   fdot+= lf*dt;
  176.   f(2:-1,2:-1)+= fdot*dt;
  177.  
  178.   return i<200;
  179. }
  180.  
  181. func laplacian(f, y,x)
  182. {
  183.   /* There are many ways to form the Laplacian as a finite difference.
  184.      This one is nice in Yorick.  */
  185.   /* Start with the two median vectors across each zone.  */
  186.   fdz= f(dif,zcen);
  187.   fzd= f(zcen,dif);
  188.   xdz= x(dif,zcen);
  189.   xzd= x(zcen,dif);
  190.   ydz= y(dif,zcen);
  191.   yzd= y(zcen,dif);
  192.  
  193.   /* Estimate the gradient at the center of each cell.  */
  194.   area= xdz*yzd - xzd*ydz;
  195.   gradfx= (fdz*yzd - fzd*ydz)/area;
  196.   gradfy= (xdz*fzd - xzd*fdz)/area;
  197.  
  198.   /* Now consider the mesh formed by the center points of the original.  */
  199.   x= x(zcen,zcen);
  200.   y= y(zcen,zcen);
  201.   xdz= x(dif,);
  202.   xzd= x(,dif);
  203.   ydz= y(dif,);
  204.   yzd= y(,dif);
  205.   area= xdz(,zcen)*yzd(zcen,) - xzd(zcen,)*ydz(,zcen);
  206.  
  207.   return ((xdz*gradfy(zcen,)-ydz*gradfx(zcen,))(,dif) +
  208.       (yzd*gradfx(,zcen)-xzd*gradfy(,zcen))(dif,)) / area;
  209. }
  210.  
  211. func pl3d(shading, z, y, x)
  212. {
  213.   /* rotate so that (zp,yp) are screen (y,x) */
  214.   /* These orientations are cunningly chosen so that the painter's
  215.      algorithm correctly draws hidden surfaces first -- see help, plf
  216.      for a description of the order cells are drawn by plf.  */
  217.   theta= 30. * pi/180.;  /* angle of viewer above drumhead */
  218.   phi= 120. * pi/180;
  219.  
  220.   ct= cos(phi);
  221.   st= sin(phi);
  222.   yp= y*ct - x*st;
  223.   xp= x*ct + y*st;
  224.  
  225.   ct= cos(theta);
  226.   st= sin(theta);
  227.   zp= z*ct - xp*st;
  228.   xp= xp*ct + z*st;
  229.  
  230.   if (!shading) {
  231.     color= [];
  232.     edges= 1;
  233.   } else {
  234.     /* compute the two median vectors for each cell */
  235.     m0x= xp(dif,zcen);
  236.     m0y= yp(dif,zcen);
  237.     m0z= zp(dif,zcen);
  238.     m1x= xp(zcen,dif);
  239.     m1y= yp(zcen,dif);
  240.     m1z= zp(zcen,dif);
  241.     /* define the normal vector to be their cross product */
  242.     nx= m0y*m1z - m0z*m1y;
  243.     ny= m0z*m1x - m0x*m1z;
  244.     nz= m0x*m1y - m0y*m1x;
  245.     n= abs(nx, ny, nz);
  246.     nx/= n;
  247.     ny/= n;
  248.     nz/= n;
  249.     color= bytscl(nx, cmin=0.0, cmax=1.0);
  250.     edges= 0;
  251.   }
  252.  
  253.   plf, color, zp, yp, edges=edges;
  254. }
  255.